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Abstract 

This paper is concerned with the identification of the geometrical structure of the 
system boundary for a two-dimensional diffusion system. The domain identification 
problem treated here is converted into an optimization problem based on a fit-to-data 
criterion and theoretical convergence results for approximate identification techniques 
are discussed. Results of numerical experiments to demonstrate the efficacy of the 
theoretical ideas are reported. 
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I. INTRODUCTION 


Domain identification problems are important in the design of engineering systems and 
frequently such problems are treated as a branch of the calculus of variations which involves 
nonlinear optimization techniques, optimal control theory, partial differential equations (el- 
liptic, parabolic, hyperbolic, etc.) and related numerical methods. Domain identification 
for elliptic systems has been studied theoretically and numerically by many authors (see 
e. g. , [5], [7], [10], [13]). For parabolic systems, a couple of numerical methods for identi- 
fying the domain or boundary have been investigated in [14], [15]. Until recently, most 
investigations concentrated on the “optimal shape design problem” which is motivated by 
numerous applications to structural, engine, airplane, ship designs, etc. (see [10] and the 
references therein). In this paper, our concern for domain identification is motivated by an 
application that is different from these shape design problems. However, as we shall see, 
the resulting theoretical aspects are closely related. Recently, associated with the use of 
fiber reinforced composite materials for aerospace structures, there is growing interest in 
the detection and characterization of large structural flaws which may not be detectable 
by visual inspection. One recent effort has focused on non-destructive evaluation methods 
(NDE) based on the measurement of thermal diffusivity in composite materials (see e.g., 
[8]). Motivated by these problems, we consider the domain identification for parabolic 
systems. 

To explain our approach, we restrict our attention to a 2-D domain identification prob- 
lem. We consider the bounded domain G(q) in two-dimensional Euclidean space as follows: 

G(q ) = {(zi,x 2 )]0 < xi < 1, 0 < x 2 < r{x l ,q)} 

where X\ — *• r(xi,q) is some parameterized real function which is assumed to characterize 
the unknown part of the boundary and q is a constant parameterization vector to be 
identified among values in a given compact admissible parameter set Q. As depicted in 
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Fig. 1.1. The spatial domain G and its boundary dG l , dG 2 , dG q , dG±. 

Fig. 1.1, we assume the boundary of G(q) consists of the following components: 
dGi = {x — (a:i, 0 : 2 ) |0 < Xi < 1, x 2 = 0} 

dG 2 = {x = {xi,x 2 )\xi = 1, 0 < x 2 < £} 

dG q = {x = (xi,x 2 )|0< xx < 1, x 2 = r(xi,?)} 

dG 4 = {x = (xi,x 2 )|xi = 0, 0 < x 2 < £} 

The measurement system is described by the following 2-D diffusion equation: 

du ^~y - dAu{t, x) + c 0 u(f, x) = 0 in T X G{q) 
dt 

with the initial and boundary conditions, 

u(0,x) = u 0 (x) on G{q) 


(1.1a) 


( 1 . 16 ) 


2 


du , , „ 

Cl dx 2 +h ( U fi = ° 

on 

T x dGi 

(l.lc) 

du 

ci \- hu = 0 

OX i 

on 

T x dG 2 

(i-i<0 

du 

dil = 

on 

T x dG q 


du , 

-ci- h hu = 0 

UX i 

on 

T x dG 4 

(ii/) 


where ci, c 0 and h are thermal diffusivity, radiation coefficient and heat transfer coefficient, 
respectively, which are given constants, and where T denotes the time interval (0, tf) during 
which the process is observed. In the above system, / is the known boundary input defined 
on T x dG\ and u 0 is the given initial function defined on Q where fl is a known bounded 
domain in R 2 such that fl D G(q) for any qeQ. The system output is assumed to be on a 
subset E of the boundary dG i, and mathematically, the observation is taken as 

y[t,x,q) = u(t,x u Q,q ) ( t,x)cT X E. (1.2) 

From a physical point of view, the system state u = u(t, x) represents the temperature 
distribution at time t at location x = (x u x 2 ) and, the boundary input / and the output 
y correspond, respectively, to the thermal source (for example, by a laser beam) and 
the observation of the temperature distribution at the surface of the material (e.g., by an 
infrared imager) (see [8] for more details). Thus, the search for structural flaws in materials 
may be formulated as an inverse problem for a heat diffusion system. The problem treated 
here is that of identifying, from input and output data {/, ub,y} on (T x dGi) x G x (T x 
E), the constant parameter vector q in Q determining the geometrical structure of the 
boundary dG q . 

In Section 2, we formulate this problem in an abstract setting in a Hilbert space. In Sec- 
tion 3, for computational purposes, we approximate the Hilbert space by finite dimensional 
subspaces and we discuss the convergence analysis for the approximate identification prob- 
lems. In Section 4, a practical optimization technique based on a finite element approach 
is outlined. Some numerical results for a simple example are given in Section 5. 
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II. PROBLEM FORMULATION AND BASIC ASSUMPTIONS 


For the discussions here, we restrict the geometrical structure of the boundary dG q by 
imposing the following hypotheses: 

(H-0) The admissible parameter set Q is a compact subset of R n \ 

(H-l) For each qeQ, we have r(g)eW^(0, 1); 

(H-2) For each qeQ, we have r(0, q) = r(l, q) = £; 

(H-3) There are constants /?i and /3 2 satisfying 0</3 1 <£</3 2 <oo such that, for qeQ, 
we have 

Pi < r i^q) < P 2 a.e. in (0,1); 

and 

(H-4) There exists a constant M such that 

l r U>ff) - r (^,9)|i,oo < M\q- q\ for q, qeQ, £e[0, 1], 

where | - |i )0 o denotes the norm of W^(0, l). 

We make the following assumptions for the class of system inputs: 

(H-5) u 0 eL 2 ( H); 

and 

(H-6) / £ £ ! (r;L ! (0,l)). 

It follows from results in ([9], Ch. 3) that, under the hypotheses (H-l)-(H-3), (H-5), and 
(H-6), for each fixed qeQ, there exists for (1.1) a unique solution u in L 2 (T-,H l (G(q))). 
Following standard procedures in optimal shape design techniques ([10], Ch. 8, p. 125), 
we introduce the affine mapping 

(z u z 2 ) = T(q)(x x ,x 2 ) 
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given by 


Zl = Xi 
_ 1X2 

Z 2 r{x u q)' 

Note that this is equivalent to xi = z u x 2 = r(zi,q)z 2 /l. Under this coordinate change, the 
system domain G(q) is transformed into the fixed domain G 

G(ff)->G = (0,l)x(0,0 

which is independent of the parameter q. Using this coordinate transformation, we obtain 
the system state u given by 

u(t,z) = u(t) o T~ X {q) =u(f,xi(zx), x 2 (z l5 z 2 )); 


this transformed state then satisfies the system equation 

E — 

.\y<2 d * 


d -^T - + «o*(M) = o in 


dzj 


i< 2 


dzj 


with 


u{0, z) = u 0 [z) oj l [q) on G 


T xG 

(2.1 a) 


(2.16) 


-7 - ■ * N f U + h(u - /) = 0 on T X dGi (2.1 c) 

r{z u q)dz 2 


du 


Cl 


dz\ 


ci .. . du 

J z * r ( 1 ^)o^ + hu = 0 


on 


T x dG 2 


(2.1 d) 


' ( 2 i ' ?) Sr “^)[ <r,( " 1, ' 7)}!+1 ^ =0 on Txaa ‘ (2 ' le) 


-ci^- + ^-2 2 r'(0, q)^— + hu = 0 on T x dG 4 . 

OZ\ t UZ2 


(2-U) 
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In this system the coefficients are given by 


a ll C 1 

(2.2) 

, V / \ c 1 r'(z u q)z 2 

an (*) = «»(*):- 

(2.3) 

a M*) := { r (z u q)y [i r '( z ^9)y4 + ^} 

(2.4) 

c i r '( 2 i »?) 

MZ) - ~ 

(2.5) 

c 1 {r'(2 1 , 9 )} 2 2 2 

Mz) - = 

(2.6) 


where r' denotes dr /dzi, while dG q has been mapped into 


dG 3 = {z = [z u z 2 )\0 < zi < 1, z 2 = £}. 


If we consider a variational formulation similar to that in [2], [9], the system dynamics can 
be described by the variational form: 


< > + CT (?) (“(<)><£) = L(t,q){<j>) for <j>zH l (G) 


(2.7) 


ti(0) = u 0 o T 1 (q) 

where the bracket < *,• > denotes the scalar product in L 2 (G ) and where cr(<7)(-,-) and 
L(q){-) denote, respectively, a sesquilinear form on x if 1 (G) and a linear functional 

on fT 1 (G). Explicitly, a(g)(-,-) and L(q)(-) are given by for <j>,xl)tH x {G) by 

^ * ] := / l l £ 2 <■«(*) J|i£ + g ^ + Si Ttjtr) dz ' 

( 2 . 8 ) 

J Q 


+ 


dz 2 , 


n ih 

L{t,q){4>) := / — rf{t,zi)[<t>\ aGl dz u 

Jo r{zi,q) 


(2.9) 
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respectively. With some tedious calculations, one can readily establish the following useful 
conditions on the sesquilinear form a. 

Theorem 1: Let \ • |y and | • |jy denote the norms in the Hilbert spaces V = i? 1 (G) and 
H = L 2 (G). Then, under the hypotheses (H-0) to (H-4), the sesquilinear forms cr(^)(-,-) 
satisfy the following inequalities: There exist positive constants /ex , A, /c 2 , and k 3 such that 
for <f>,il)eV we have 

<7(9) (^>$) > kil^lv ~ (2-10) 

<r{q){<t>,'l>) < fcaMvMv (2-11) 

- *(5) (^.0)1 < k 3 \q- q\\(f>\v\^\v for all q,qeQ. (2.12) 


Proof: We wish to show first that the sesquilinear form a is coercive. From (2.2)-(2.4) 
and (2.8) , the principal part of the differential operator becomes 


£ <>«(*, «)66 = c . 6 * - 


i,j< 2 


r{zi,q) ^ {r(zi,q)V 

for (6,6)ei? 2 . 

By simple calculations and from (H~3), we have 

Cl* l 161 s . 161 s 


(2.13) 


y~! a i,j[z,q)£itj > 


«,J<2 


f lii 

U r '( 2 i>9)) 


+ 


9 )} 2 z! + £2 {^(21,9)} 


)* 


^i(|ei| 2 +l6| 2 ) (2.14) 


where 




Ci^ 2 

2/?f * 


This means the operator is strongly elliptic. For the coefficients bj(z,q)(j = 1,2), from 
(2.5) and (2.6), it follows that 


\bi{z, 9)| < j- sup |r'(*i,9)| 

Pi zido.il 

(2.15) 

|6 2 (^,9)l < ^ sup |r'(2 1)9 )| 2 

Pi *l«[0,l] 

(2.16) 
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respectively. We note that, from (H-0) to (H-4), 

sup |r'(a,?)| < R (2.17) 

where R is some constant independent of q. Hence, we obtain 


\bj{z,q)\ < K 2 < oo for j — 1,2, 


(2.18) 


where (assuming R > 1) 


K ,= 


ci LR 2 

~1T 


(2.19) 


For the last two boundary integrals in (2.8), by virtue of (H-3), the following inequality 


holds: 

r 1 PI j rl r 

I r[^ 2 ]aci^i + / h[4> 2 ] aG3 udGidz2 > h / _[(f> 2 ] aG ds (2.20) 

JO J 0 J dG 

where ds denotes a line element on dG and dG = dG\ U dG 2 U dG±. From (2.14), (2.18), 
and (2.20), we can derive the coercivity property of the sesquilinear form. Namely, the 
sesquilinear form satisfies 


>Ki j “ K 2 




( 2 . 21 ) 


where 


K > = min [w h , 


Friedrich’s second inequality ([1], p. 124) asserts that if dG is a nontrivial subset of dG, 
then there exists a positive constant a such that 



(2.22) 
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By applying this to the last parenthesis of (2,21), we can conclude that 




(2.23) 

where 

K-Z + KHCfU 

4 K i 4 


respectively. 



To prove 

the boundedness of o(q), we note that 



k(?)(<M) 1 <12/ ja ai ^ z,q ^dzd Zi dz ^ 
*,y<2 3 

+12 / f d b }( z >Q)j^-y dz \+ 1 / f 6 coH>dz\ 

rl lh 

+| y o + y o M00]aG 9 uaG 4 <fea|- 

(2.24) 

The first three integrals of RHS in (2.24) satisfy 



fjg M z *v)\ \4W\t\v 

,>J - 2 J ’ * £ (0,’llx[0,t) 

(2.25) 


fG bj ( Z,q ^^z : ^ dZ \ ~ 2 l 6 i'( Z ’9)M^|v|0k 

J - 2 3 «[ 0 ,lTxl 0 ,t] 

(2.26) 


1 J c 0 <t>xl)dz\ < c 0 \(j>\v\ip\v. 

(2.27) 

respectively. 

From (2.2) to (2.4), it follows that, under hypotheses (H-l) and (H-3) (see 

( 2 - 17 )). 

sup \a,ij(z,q)\ < Kq < oo 

(2.28) 


**y< 2 

*c(0,l]X[0,t] 


where 




Ke=^-{R 2 + l). 

From (2.18) and (2.25)-(2.28), we have 

rl rt 

\o{q){<j>,i>)\ < (4 Kq + 2 K 2 + Co) \<f>\ v \4>\ v + I y o — -^[^i^aG^i + J q h[</>^]aGr 2 uaG 4 ^ 2: 2 | 

(2.29) 
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Furthermore, the boundary integral term satisfies 



th 

r{zi,q) 





< K 7 \4>\v\fp\v 


which follows from the trace inequality ([l], p. 124) 


L 


-W ]aa ds < 7{ G )\<t>\ 

OG 


2 

V) 


(2.30) 


(2.31) 


where 


K 7 = 


hii{G) 

(h 


Consequently, we can prove the boundedness property of a{q) (•, ■). 

To establish the continuity property, we note that, for any q and qeQ, 

\o{q){<i>,^) - °{q){<t>,'I>)\ 

+1 Jf- £(M?) - *>;(9))||-^l 

J J G j< 2 oz j 

//glMs) -“»(«)! l|^l 1^1* 


+ 


+ 


f j 6 \bi{q) ~ h(q)\ \~\ \ip\ dz + J J & \b 2 {q) - b 2 (q)\ \^-\ \ip\di 


-(<?) !■(?) 

Under the hypotheses (H-l) and (H-3), we argue that 

c\t R 

|ai 2 (g) - o 12 (g)| < , 1) |r(g) - r(?)|i >0O 

|a 22 ( 9 ) - a 22 {q) I < — l)l r (?) - r(g)|i,c 

Pi Pi 

IM?) ~bi(q)\ < ^-max{^-,l)\r{q) -r(9)|i,oo 
Pi Pi 


(2.32) 
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2cilR_ _ ( RPi 


M?) - *2(5)1 < l)l r (?) “ r(5)|i,o 


and 


| ^rKi) | -« ||r(?) ' r ® ||o( °’ , >- 

Applying these inequalities into (2.32), we have 

\°{q){4>A 0 - I < Kg\r(q) - r($)|i |OO |^|y|0|y 


(2.33) 


where 

__ 2ci£ .R . 2c\t? R .2R02 . c\ . R . 2c\LR f R0 2 

if, = — ma I (^,l)+^ ? -m aI (-^-,X) + -m«x(- 1 l)+-g r max(-^-,l) + ^ 7 (G). 

From the hypotheses (H-4), we can thus infer the continuity of the sesquilinear form 
o(g)(-, •) with respect to the parameter q in Q. The proof has been completed. 

For the system (2.7), the output can be represented as the restriction of u(t) to a subset 
E C dGi of positive measure, i.e, 


y{t,q) = u(f,g)| E . (2.34) 

We assume (see (H-0)) throughout that the admissible parameter set Q is a given compact 
subset of R n . The fundamental identification problem considered here is based on the 
fit-to-data functional (see [2]) given by 

J W = \J 0 ' 9) - Vd{t) || \dt (2.35) 

where F — L 2 (E), {yd{t)}uT are given observed data, and y(t,q) is the solution of (2.7) 
corresponding to qeQ. Then our problem is stated as follows: 

(IDP) Find q*eQ which minimizes J(q) given in (2.35) subject to the system (2.7) and 
(2.34). 

In the next section, we consider a family of approximating identification problems associ- 
ated with (IDP). 
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III. APPROXIMATE IDENTIFICATION PROBLEMS 


The approximation scheme we have employed is based on the use of a finite element 
Galerkin approach to construct a sequence of finite dimensional approximating identifica- 
tion problems. Let us choose Ujv=i{^f^}£Li as a set of basis functions in H 1 (C). That is, 
for all N, a re linearly independent and Un span^^}^ is dense in the V norm 

in V = J9 rl (G). We choose the approximation subspaces as 

H n := span{4 >^ , <f>% , ■ ■ • , 

Then, we can define the approximate solution of Eq. (2.7) by 

u N (t,q) = (3-1) 

i<N 

where i v^(t,q) are chosen such that for j — 1,2, • • • , N, 

< (3.2a) 

and 

« N (0) = 53 < “o ° T -1 (?),<?!)f > cj>i . (3.2 6) 

i<N 

Hence the system (2.7) and the output (2.34) can be approximated by solving the system 


where 


C N w N (t,q) + A N (q)w N (t,q) = F N (t,q) 
w N (0) = JUq 

y N (t,q) = 53 (*>g)[^fb 

i<N 

[C N h =< < > for t\ j = 1 , 2 , • • • , N 

[^(g)].,; := o(q)(<f>f,<f>^) for *,J = 1,2, • • • ,N 

[w N {t,q)]i := w?(t,q) for t = l,2,---,JV 
[**(*,*)], := for i = 1,2, ■ ■ ■ , JV 
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(3.3a) 

(3.36) 

(3.4) 



Wo]i :=< «o o r l {q),<l>? > 


for * = 1,2,- • • , N. 


The approximating identification problems thus take the following form: 

(AIDP) N Find q N eQ which minimizes 

j"( t ) = \f o 'W(t, q) -ymUt (3.5) 

subject to the approximating system (3.3) and (3.4). 

Our convergence results for the finite element schemes are summarized in the following 
two theorems. 

Theorem 2: Let {q M } C Q be a sequence such that q M — ► qeQ as M — *• oo and 
let u N (q M ) and u(q) be the solutions of Eqs. (3.3) and (2.7) corresponding to q M and q, 
respectively. Then, under hypotheses (H-0) to (H-6), we have u N (q M ) —> u{q ) strongly in 
L 2 {T]H 1 {G)) as N,M-+o o. 

Theorem 3: Let q N be a solution of the problem {AIDP) N . Then the sequence {q N } 
admits a convergent subsequence {q Nk } with q Nk — ► q as k — ► oo. Moreover, q is a solution 
of the problem (IDP). 

The proof of Theorem 2 follows from the general convergence framework for parameter 
identification problems given in [3] and [4]. To ensure the desired convergence, it suffices to 
show that the sesquilinear form a(^)(-,') satisfies the continuity, coercivity and bounded- 
ness conditions as stated in [3], [4]. But this is a result of Theorem 1 under the hypotheses 
(H-0) to (H-4) . 

The proof of Theorem 3 can be carried out by using Theorem 2 and the compactness 
of Q. Since q N is a solution of the problem ( AIDP) N , it is clear that 

J N (q N )<J N (q) for WqeQ. (3.6) 

Thus, if we can argue that for any q M — > q in Q, 

y N {q M )-+y(q) in L 2 {T]F) as ]V,M^oo, 
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then, we can obtain the desired inequality 


J{q) < J{q) f° r VgeQ 

by taking limits in (3.6). But the needed arguments follow immediately from Theorem 2 
since 

IliW) - lf(4)|ll W ) < if||6 K (4''*) - S(4)lli.(T!V) 

where K is independent of q Nk and q. 

IV. OPTIMIZATION TECHNIQUES FOR THE APPROXIMATE ESTIMA- 
TION PROBLEMS 

Let q N be an optimal solution of the problem (AIDP) N . Then a necessary condition 
for q N to be optimal is characterized by 

V g J N (q N )-(q-q N )>0 for VqeQ (4.1) 

where V, denotes the gradient of J N (q) with respect to q. From Eq. (3.5), we have for 
k = 1,2, • • • ,n 

[V,J W (?)1» = /"K(t,4))'(C,V(<,,) - y/W)* 

J o 

where 

=< <}>?, <t>? >f for », j = 1, 2, • • • , N 

[y/(t)L- =< >f for j = l,2,--.,iV. 

Using the same procedure as in [9], we can evaluate the gradient vector by (for k = 

1,2, • ■ ■ ,n) 

(V, •/"(«)]* = f' v"(t, )]«-"((,?) - V n F»(t, g )}dt (4.2) 

Jo 

where v N (t,q) is the solution of the adjoint equation, 

- C N v N (t , q) + A* N {q)v N (t, q) = Y d N {t) - C*w N {t, q) (4.3a) 
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»"(*/, ?) = 0. 


(4.36) 


In Eq. (4.3a), the matrix A* N (q) is given by 

where G*(q){-,-) is the adjoint sesquilinear form of <j(g)(’,-) defined by 


--IL 


a* {q)(<f>, 0) 

E <>«(«.*)## - +(<=.- E §^(«.*)W 

i,j<2 ° z i OZ ' J<2 J<2 47 Z J 


dz 


/•l £/! rl 

+ J o r( g ^ WAi + J o K^dGiudGidzi 

'Chsr'^S, 

Consequently, the optimality condition (4.1) of the problem ( AIDP) N can be characterized 
by 


E ?'')'{[V„A''(5'')]u. A '((,j'')-V !l F"((,5' V )>(9*~9 4 '') >0 for all ,£«. 

(4.4) 

In the sequel, we discuss computer implementation of numerical schemes for the prob- 
lem ( AIDP ) N . Since we can evaluate the gradient of the cost function using (4.2), many 
optimization techniques for the constrained problems are readily applicable to our problem 
(see [11] and the references therein). For ease in exposition, here the compact set Q C R N 
is assumed to be defined by 


Q = {q = (q l ,q 2 ,‘--,q n )^R N \^q < q} (4.5) 

where II and q denote a given constraint matrix and vector, respectively. For the numer- 
ical results reported in this paper, we used the gradient projection method [12] which is 
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a particularly useful technique for optimization problems with the linear inequality con- 
straints such as those given in (4.5). We use this method as presented in [12]; the iterative 
algorithm for finding q N can thus be stated as follows: 


Step 0: Choose an initial value q^ in Q and set i = 0. 

Step 1: If IlgW < (j se t 

g (i) = -V g J N (qM) 

and proceed to Step 8; otherwise, proceed to Step 2. 

Step 2: Compute the current direction by 

M = fV . JJ, (9 (i) ) 

5 |PV,J»(«(f))| 

where 

p = / - n;(n,n;)- l n, 

and n p includes the gradient of all currently active constraints associated with matrix IT. 
If g ( l ) ^ o, proceed to Step 5; otherwise, proceed to Step 4- 

Step 8: Compute A^ in satisfying 

J N (q (<) + A % in 9 {i) ) = min J N {fl + A </«) 

A £ [o,A] 

where A is the largest step that may be taken from q M along gW without violating any 
constraint. If A^ tn = A, then add the new contraints to the matrix n p and proceed to Step 
4’, otherwise, the new approximation to the solution is given by 

,«■»> = 5 » + aJIjW. 

Replace t + 1 by i and return to Step 1. 

Step 4-' Compute the vector 6(q) by 

%») = -(n p ny-'n p v p j(,W). 
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If all components of 6 are nonnegative, then set 


and terminate the computation; otherwise, delete the column of n p corresponding to the 
smallest component of 0(<jM), replace * + 1 by i and return to Step 1. 


V. NUMERICAL PROCEDURES 


In a series of numerical experiments, we used a test example constructed as follows: 
We chose a function r(g), generated the corresponding solution numerically, added random 
noise, and then used this as “data” for our inverse algorithm. The parameter function 
r(£,q) to be identified is a piecewise cubic polynomial function (see [6] for more details). 
We denote the knot sequence for r by 


and the unknown function r(£,g) is given by 

r (£>0 = P?(0 

= a l,i + a 1,i{£ ~ r r) + a 3,»'(£ ~ r t”) 2 /2 + a 4,t(£ — t «") 3 /6 (5.1) 

for r" < £ < x *' = 0,1, 

The unknown parameter vector q = {&}£_! is then given by 

= r(r,) for * = 1,2, •••,«. (5.2) 


Further, we assume 

Po(0,q) =Pn(l,9) = t (5.3) 

Po(0,9) =p' n (hq) = 0. (5.4) 

Substituting (5.2), (5.3), and (5.4) into (5.1), the coefficients {a*,} can be determined 
uniquely and r(£, q) satisfies the hypotheses (H-l), (H-2), and (H-4). In order to guarantee 
the hypothesis (H-3), we impose the constraints 


Pi<Qi< Pi for t = 1, 2, • • • , n. 
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Hence, the matrix n (2 n x n) and the vector q (2 n x 1) defining the admissible parameter 
class Q (see (4.5)) is given by 


1 

-1 


n = 


1 

-1 


0 




A 

-A 


A 

-A 


To discretize the system model by the finite element method, the domain G is divided 
into a finite number of elements {ek}Jf=i ( K < N ) and a number of nodes defined by 
{zi = (z[, are selected in G. For convenience of computations, we set t = 1 in G. 

Each element is preassigned as an axiparallel rectangle with nodes at the vertices. The 
restriction of <f>^ to any element e* is given by the bilinear polynomial form, 


= 4* + + C S Z 2 + ^ziz 2 (5.5) 

for z = (zi,Z2)eek A: = 1,2,*--, JK" and » = 1,2, • • • , N. 


The coefficients {cj^} can be chosen such that each polynomial form (5.5) satisfies the 
properties of a piecewise bilinear basis function (see e.g., [l], Ch. 5). The integration of 
element matrices C N , A N (q), A*^?), and Cjf, and the element vectors F N and Yj* can 
be computed numerically by a Gauss-Legendre formula. Thus, the state model (3.3) and 
its adjoint system (4.3) can be solved numerically by an implicit scheme with respect to 
discrete time t = ih (t = 0, 1, • • ■ , m) where h — tf/m. The evaluation of cost functional J N 
and its gradient W q J N is the computationally expensive part of our algorithm since these 
involve the integration of the states w N (t,q) and the adjoint states v N (t,q) with respect 
to time t over T. This can be accomplished by using the two-point Gauss formula. 

The input data are preassigned as 


u 0 {z) -- -10, 


for zeG 


/(0 = 0, for e^dGi. 
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The known parameters c l5 c 0 , and h in Eq. (1.1) were set as 


Ci = 0.034, c 2 = 0.001, h = 0.1. 

The observed data {y<z(t)} were generated by solving the finite element model (3.3). The 
number of finite elements and nodes in the numerical experiments were set as K = 256 (= 
16 x 16) and N = 289 (= 17 x 17), respectively. The final time and number of time 
divisions were taken as tj = 10 and m = 100. Random noise at various levels from 0% to 
50% was added to the numerical solution, thereby producing simulated noisy “data” for 
the algorithm. The set E relative to data acquisition was given by 

S = U w., 

r- l 

where N Xj denotes a neighborhood of points Xj at dG i, i.e., 

N x . = ( X j-e,Xj + £) for j= 1,2, -,p. 

Using such data, the estimation algorithm given in Section 4 was tested. 

Example 1: In this example, the dimension of unknown vector was taken as n = 4 and the 
knot sequence was given by 

T* = i/ 5 for * = 1,2, •••,5. 

The values of the true parameters were chosen as 

qi = r(r*) = 0.8 for t = 1,2, 3, 4. 

The lower and upper bounds of the unknown parameter vector were taken as /?! = 0.3 and 
/?2 = 1.1, respectively. The initial guesses for the parameters were given by 

qj 0) = 1 for t = 1,2, 3,4. 

The number of sensors was taken as p = 9. Table 1 shows the estimated parameter numer- 
ical results for the data with noise free, 5%, 10%, and 50% relative noise and Figure 5.1 
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shows the estimated parameter function r(£, q N ) and true function r(f , q) which correspond 
to the estimated boundary shape and true boundary for the 10% noise case. 



A 

9i 

A 

92 

A 

93 

A 

94 


9. ~ 9» 

3 

True Value 

0.800 

0.800 

0.800 

0.800 



Initial Guess 

1.000 

1.000 

1.000 

1.000 



Noise 

iteration 6 

0.850 

0.882 

0.880 

0.849 


Free 

iteration 13 

0.820 

0.819 

0.821 

0.819 



iteration 17 

0.810 

0.785 

0.810 

0.806 



iteration 6 

0.851 

0.879 

0.878 

0.844 

■ 


iteration 13 

0.828 

0.821 

0.818 

0.820 

1 

m 

iteration 17 

0.815 

0.797 

0.791 

0.814 

1 

10% 

iteration 6 

0.849 

0.853 

0.861 

0.834 

HU 

Noise 

iteration 14 

0.787 

0.847 

0.835 

0.750 



iteration 18 

0.827 

0.792 

0.799 

0.796 

. 

50% 

iteration 7 

0.922 

0.820 

0.748 

0.808 

3.36 X 10~ 2 

Noise 

iteration 15 

0.902 

0.793 

0.772 

0.886 

3.40 x 10- 2 


iteration 19 

0.813 

0.783 

0.727 

0.794 

1.90 X 10“ 2 


Table 5.1. True Value and Estimated Values in Example 1. 


Example 2: We chose the same dimension of unknown parameter vector as in Example 1 
and we also used the same knot sequence. In this example, however, the values of the true 
parameters were preassigned as 

9i = 94 = 0.9 


and 


92 = 93 = 0 . 6 , 


respectively. The lower and upper bounds, initial guess of unknown vector, and number 
of sensors were given by the same values as in Example 1. Table 5.2 shows the numerical 
results obtained here for the various sets of noisy data. Figures 5.2 and 5.3 represent the 
estimated parameter function for the case of 20% and 50% noisy observation case. 
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Figure 5.1. True Function and Estimated Function in Example 1 (10% Noise). 


21 



A 

9i 

92 

93 

A 

94 


H 

3 

True Value 

0.900 

0.600 

0.600 

0.900 



Initial Guess 

1.000 

1.000 

1.000 

1.000 



Noise 

iteration 16 

0.817 

0.792 

0.778 

0.821 


Free 

iteration 25 

0.894 

0.607 

0.602 

0.893 


5% 

iteration 16 

0.953 

0.694 

0.811 

0.906 


Noise 

iteration 26 

0.896 

0.604 

0.605 

0.898 


10% 

iteration 16 

0.902 

0.807 

0.674 

0.949 


Noise 

iteration 27 

0.908 

0.594 

0.581 

0.896 


20% 

iteration 16 

0.917 

0.683 

0.812 

0.904 


Noise 

iteration 27 

0.902 

0.603 

0.610 

0.887 


25% 

iteration 16 

0.939 

0.819 

0.684 

0.962 

■ 

Noise 

iteration 26 

0.868 

0.563 

0.563 

0.874 

mEBBSM 

50% 

iteration 16 

1.02 

0.618 

0.680 

0.915 


Noise 

iteration 28 

0.951 

0.574 

0.599 

0.953 



Table 5.2. True Value and Estimated Values in Example 2. 


Example 3: In this example, we deal with a somewhat more difficult case as compared 
with Examples 1 and 2. We set the dimension of parameter space as n = 8 and we chose 
the knot sequence as 

{r-}-=o r* = i/9 for * = 0, 1,2, • • • , 9. 

True parameter values were given by 


9i = 9s = 0.99, 


92 = 97 = 0.98, 


93 = 94 = 0.94, 


and 


95 = 96 = 0.60, 


respectively. Figure 5.4 shows the corresponding boundary shape to be identified. The 
number of sensors was taken as p = 17. The bounds and initial guesses for the parameter 
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Figure 5.4. Unknown boundary shape in Example 3. 

vector were the same as in Examples 1 and 2. We ran numerical experiments for the case 
of noise free, 5%, 10%, 20%, and 50% noisy observations. Table 5.3 shows the estimated 
parameter vector obtained here. Figure 5.5 represents the estimated boundary curve for 
the 10% noisy data. 
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A 

9x 

A 

92 

A 

93 

A 

94 

A 

95 

A 

96 

A 

97 

A 

98 


na 

■ 

True Value 

0.990 

0.980 

0.940 

0.600 

0.600 

0.940 

0.980 

0.990 



Initial Guess 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 



Noise 

iteration 16 

1.039 

1.027 

0.887 

0.808 

0.674 

0.926 

1.029 

1.033 

SI 

Free 

iteration 23 

1.014 

1.007 

0.943 

0.624 

0.615 

0.943 

1.003 

1.017 



iteration 16 

1.042 

1.031 

0.879 

0.749 

0.666 

0.920 

1.049 

1.031 

1 


iteration 23 

1.023 

0.987 

0.948 

0.626 

0.628 

0.937 

1.002 

1.020 

1 

10% 

iteration 16 

1.037 

1.058 

0.899 

0.733 

0.721 

0.881 

1.044 

1.025 


Noise 

iteration 24 

1.010 

1.009 

0.955 

0.586 

0.578 

0.948 

0.994 

1.000 


20% 

iteration 16 

1.034 

0.989 

0.950 

0.592 

0.595 

0.994 

1.018 

1.061 


Noise 

iteration 24 

1.022 

0.985 

0.915 

0.616 

0.610 

0.941 

0.993 

1.061 


50% 

iteration 16 

1.090 

1.176 

0.872 

0.637 

0.654 

0.933 

1.166 

1.096 


Noise 

iteration 24 

1.030 

1.033 

0.939 

0.566 

0.603 

0.933 

1.046 

1.052 



Table 5.3. True Value and Estimated Values in Example 3. 


Throughout the numerical experiments, we checked the robustness of the algorithm with 
respect to noise in the observed data. Results in three examples indicated that the algo- 
rithm worked very well (i.e., as expected) for various noise levels. Furthermore, we checked 
the sensitivity of the algorithm with respect to the number of sensors. Specifically, we com- 
pared in Examples 2 and 3 the number of sensors (p) with the dimension of parameter 
space (n). In Example 2, for data with p = 5(>n = 4), the algorithm still yields an almost 
identical fit (to that for p = 9) even in 50% noise case while the fit could not be achieved 
under the reduced observation case p = 3(< n). Also, in Example 3, (where n = 8) the 
fit could not be obtained with p = 3 or p = 5, while the algorithm performed well with 
p = 9 (> n). Carrying out a large number of other numerical tests in addition to those 
reported for Examples 2 and 3, we suggest that the algorithm requires a number of sensors 
which is at least equal to the number of dimensions of parameter space, i.e., p > n. 


VI. CONCLUDING REMARKS 


In this paper, we have discussed techniques for estimating the system boundary shape 
in two dimensional parabolic systems. By using a simple coordinate transformation tech- 
nique, the parabolic PDE defined on unknown spatially varying domain was converted 
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Figure 5.5. True Function and Estimated Function in Example 3 (10% Noise). 
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into the same type PDE with unknown coefficients defined on a fixed domain. Thus, our 
fundamental approach was placed within the theoretical framework for parameter identifi- 
cation problems given in [2], [3], and [4]. The practical utility of our algorithm is supported 
through a series of numerical experiments, a summary of which is given in Section 5. These 
simulations were carried out on the Sun Microsystems at ICASE, NASA Langley Research 
Center. For three different numerical examples, using data with no noise, the proposed 
algorithm yields an almost perfect fit, while, as expected, the fit degenerates significantly 
as noise in the observation becomes more pronounced. 

Although here we discuss only the case where the unknown boundary shape is rep- 
resented by a simple function of one variable, our basic parameter estimation ideas and 
techniques can be readily extended to consider more general classes of geometrical struc- 
tures for the system boundary. For example, we may also treat the case where the unknown 
boundary shape is characterized by 

r{q,xi,Xi) = 0 for [x u x 2 )eR 2 . 

We are currently pursuing investigations for these cases. 
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